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RESUMEN 


Se planteó la necesidad de estudiar los conceptos fundamentales de la simulación mediante el 
desarrollo de un algoritmo basado en diferencias finitas, para simular la evolución espacial y 
temporal, de un fenómeno físico particular (conducción del calor) modelado con ecuación 
diferencial parcial parabólica. Este artículo es el primero de una serie de dos, en el cual se 
verifica la validez del algoritmo diseñado para 1D, mediante la utilización de resultados 
obtenidos en un artículo base como patrón de comparación; en el segundo artículo se presentará 
el estudio del problema en 2D tomando como base en algoritmo y la metodología planteada 
aquí. 
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ABSTRACT 


The necessity considered to study the fundamental concepts of the simulation by means of the 
development of an algorithm based on finite differences, to simulate the space and temporary 
evolution, of particular a physical phenomenon (conduction of the heat) modeled with equation 
parabolic partial differential. This article is first of a series of two, in which the validity of the 
algorithm designed for 1D is verified, by means of the use of results obtained in an article it 
bases like comparison pattern; in the second article the study of the problem in 2D will appear 
taking as it bases on algorithm and the methodology raised here. 


Keywords: Simulation, algorithm, finite differences, space and temporary evolution. 
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1. INTRODUCCIÓN 


Se planteó la necesidad de desarrollar un algoritmo basado en diferencias finitas, para 
simular la evolución espacial y temporal, de un fenómeno físico particular (conducción del 
calor) modelado con ecuación diferencial parcial parabólica, que permita estudiar los conceptos 
fundamentales de la simulación. Este artículo es el primero de una serie de dos, en el cual se 
verifica la validez del algoritmo diseñado para 1D, mediante la utilización de resultados 
obtenidos en (Millán et al., 2011) como patrón de comparación; en el segundo artículo se 
presentará el estudio del problema en 2D. 


I. ANTECEDENTES DEL TEMA. 


Se toma como base el artículo (Millán et al., 2011) en el cual se estudia el problema de 
la difusión de calor en 1D, se analizan los problemas asociados a un método explícito para la 
solución numérica en diferencias finitas de la ecuación diferencial parcial parabólica. Se 
obtienen resultados que corroboran el límite de estabilidad teórico y se verifica la convergencia 
mediante la comparación con resultados proveniente de la solución analítica, se manipulan la 
ventana de tiempo de simulación, el paso espacial y temporal. No se da ningún detalle del 
software utilizado refiriéndose a el sólo como “científico”. 


El trabajo de investigación presente se relaciona con el antecedente (Millán et al., 2011) 
en que desarrollamos un algoritmo para simular el mismo fenómeno y lo validaremos con los 
resultados del antecedente, los resultados obtenidos y el análisis de los fundamentos de la 
simulación presentados en el presente trabajo, sirven de base para un segundo artículo en donde 
se abordará el problema en 2D. 


II. OBJETIVO DEL TRABAJO. 


El objetivo de este trabajo es estudiar los conceptos fundamentales de la simulación 
numérica, mediante algoritmo basado en método explícito de diferencias finitas, que desarrolle 
la evolución espacial y temporal del fenómeno de la conducción del calor modelado con 
ecuación diferencial parcial parabólica. 


III. FUNDAMENTO DE ESTUDIO. 


El objeto de estudio del presente trabajo es el fundamento de la simulación numérica, 
específicamente el refinamiento de mallas, condiciones iniciales y de bordes, convergencia y 
estabilidad del algoritmo de solución de la ecuación diferencial parcial aproximada mediante 
diferencias finitas y método explícito. Se utiliza como fenómeno físico a estudiar la conducción 
de calor en un material homogéneo e isótropo con una fuente de calor exógena al material que 
proporciona el patrón de calentamiento inicial. Tabla 1 
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TABLA I 


PROPIEDADES DEL MATERIAL OBJETO DE ESTUDIO 


Conductividad Calor Difusión 
Densidad 
térmica específico térmica 
Pp k E a 
7.8 g/cm? 0.13 calV/s.cm.*C 0.11 calV/g..C 0.15 cm?/s 
IV.METODOLOGÍA. 


Se ejecutaron nueve pasos orientados por los siguientes objetivos específicos: 


1. Estudiar el fenómeno físico a simular y el material utilizado en (Millán et al., 2011). 

2. Plantear la ecuación en diferencias finitas para el método explícito y el patrón inicial de 
temperatura incluyendo los bordes (Millán et al., 2011). 

3. Establecer el núcleo del algoritmo iterativo a desarrollar en este trabajo. 

4. Establecer el límite teórico de estabilidad CFL según las características del fenómeno y del 
material estudiado. 

5. Diseñar el algoritmo e implementación en Matlab” 

6. Validar el algoritmo utilizando como patrón datos reportados en (Millán et al., 2011). 

7. Afinar la ventana temporal de simulación (T) 

8. Afinar Ax. 

9. Afinar At. 

V.DESARROLLO 

MÉTODOS Y MATERIALES. 
El trabajo se desarrolló en nueve pasos, cada uno orientado por un objetivo específico: 

1. Estudiar el fenómeno físico a simular y el material utilizado 
La conducción del calor en una dimensión espacial en un material homogéneo e isótropo, 
que no intercambia calor con el ambiente se modela matemáticamente mediante la 
siguiente ecuación diferencial parcial (Romero et al., 2001): 

ou 0*tu k 
Ot ox? pc 

Donde a. depende del material y es denominado coeficiente de difusión del calor, k es el 
coeficiente de conductividad térmica, pes la densidad y c el calor específico (pp. 95-98 ). 
Se tomó como material para el estudio el acero descrito en (Millán et al., 2011). 

2. Plantear la ecuación en diferencias finitas para el método explícito y el patrón inicial de 


temperatura incluyendo los bordes (Millán et al., 2011) 
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Se aproximó la derivada temporal con una diferencia progresiva dada por la ecuación 
(Romero et al., 2001). 


Ou _ u(x,t + At)- u(x,t) 


2 
Ot At Se 


Para la aproximación de la derivada espacial se utilizó la diferencia de punto central, para 
en método explicito está dada por (3). 


9? + Ax,t)- 2u(x,t) + =Ax,t 
pe La u(x E) Aia u(x x,t) 3) 
OX Ax 


El patrón inicial de temperatura incluyendo los bordes, está dado por las ecuaciones (4) y 


(5) 


u(X,0)= 200-10X  10<X<20 (5) 


3. Establecer el núcleo del algoritmo iterativo. 
El núcleo del algoritmo iterativo se obtiene vinculando t con j y x con i. 


Para el método explicito se igualan (2) = (3) y se reescriben utilizando notación matricial dando 
como resultado la ecuación (6). 


u(j+L0)= Au i+D+ (1-22) 0+2uji-D), Zo 


4. Establecer el límite teórico de estabilidad CFL según las características del fenómeno y del 
material estudiado. 


Si observamos la ecuación (6), el coeficiente A = a(At/Ax?) es determinante en el avance 
del algoritmo hacia el estado estable. En (Ezquerro,2013)se establece que para garantizar la 
estabilidad el coeficiente debe cumplir con 0 < A < 0.5, lo cual implica que el límite de 
estabilidad se da cuando At/Ax* = 1/(2a. ) (p. 98). 


La simulación se iniciará con un paso Ax =2 y At = 0.5 lo cual implica un 1=0.018 esto 
predice un comportamiento sin que aparezca el fenómeno de la inestabilidad. 


5. Diseñar el algoritmo e implementación en Matlab”. 
El núcleo del algoritmo iterativo es la ecuación (6), y será implementado en el lenguaje de 


programación de Matlab”. MATLAB Version 7.4.0.287 (R2007a) 
Operating System: Microsoft Windows XP Version 5.1 (Build 2600: Service Pack 3) 
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6. Validar el algoritmo. 


Se realizó la validación mediante el contraste con los resultados obtenidos en (Millán et al., 
2011) con la ventana temporal T=120 para Ax=2 y At= 0.5. 


7. Afinar la ventana temporal de simulación (T). 


Estudiar el comportamiento del algoritmo al afinar la ventana temporal de simulación (T) 
dejando fijos At y At, midiendo el acumulado L1. 


8. Afinar Ax. 


Estudiar el comportamiento del algoritmo al afinar Ax dejando fijos At y T. 


9. Afinar At. 


Estudiar el comportamiento de la evolución temporal al afinar At dejando fijos Ax y T. 
VIPRESENTACIÓN Y DISCUSIÓN DE RESULTADOS. 


1. Algoritmo implementado en Matlab”. 

% Define dominio e inicializa la Grilla 
L=20;  %kExtensión del dominio en x (cm) 
T=125;  —% Extensión del tiempo en (s) 
Dx=2;  %.Discretización en Xx 

x= 0:Dx:L; % Coordenada x 

Dt= 14; % Discretización en t 

t= 0:Dt:T; % Coordenada t 

% Grilla 


Nx = round((L/Dx)); % Número de nodos en x 

Nt = round((T/Dt)); % Número de nodos en t 
i=1; % Iteración en la componente espacial x 
j=1; % Iteración en el tiempo 

u(1,1)=0; % Inicializa la Grilla 


forj= 1:Nt+1 
fori= 1:Nx+2 
u(j,1)=0; % Limpia la Grilla 
end 
end 


% Patrón inicial de temperaturas 
ui(1,1)=0; 
for i=1:Nx+2 
if (i-1)*Dx <= 10 
ui(1,1)=10*(i-1)*Dx; 
end 
if (i-1)*Dx > 10 
ui(1,1)=200-10*(i-1)*Dx; 
end 
end 
u(1,:)=ui; 


%:Gráfica Condiciones iniciales 
vig(1,1)=0; 
for i=1:Nx+1 
uig(1,i)= ui(1,i); 
end 


Ciencia e Ingeniería 2014 Vol. (1) N* (1) 


Lyon E. y Rosales L. 


FUNDAMENTOS DE SIMULACIÓN MEDIANTE ALGORITMO EN DIFERENCIAS FINITAS 
DE LA EVOLUCIÓN ESPACIAL Y TEMPORAL EN LA CONDUCCIÓN DEL CALOR 


plot(x,uig(1,:),'1) 


% Característica del material 

c=0.11; % Capacidad térmica cal/(g.*C) 

p=7.8; %Densidad g/(cm3) 

k=0.13; % Conductividad térmica cal/(s.cm.*C) 
a=k/(c*p); % cm2/s 

Y =a*Dt/(Dx*Dx);% Factor de convergencia y estabilidad 


% Simulación de patrón de temperaturas 
forj=1:Nt 
fori = 2:Nx 
UG +Li)=Y*(uG Du i-D) 1-2); 
end 
end 


% Temperaturas a graficar 

ug(1,1)=0; 

forj= 1:Nt+1 
fori=1:Nx+1 
ug(j,i)=u(j,1) ; % Grilla con puntos a graficar 
end 

end 


%Gráfica de la distribución en estado estable 
hold on 


plot(x,ug(Nt,:),'g') 


%Gráfica de la evolución temporal 

t= 1:12:120; 

fori= 1:Nx 
plot(t,ug(:,i),'1'); % Grilla con puntos a graficar 
hold on 

end 


2. Validación del algoritmo. 


Se realizó la validación mediante el contraste entre los resultados de la solución numérica 
obtenida en este artículo y los resultados obtenidos en (Millán et al., 2011) con la ventana 
temporal T=120 para Ax=2 y At= 0.5. 

A continuación se presentan las gráficas de: las condiciones iniciales, la solución numérica 
y la solución analítica reportadas en (Millán et al., 2011) 

En la GRÁFICA 1 puede observarse que para las condiciones impuestas de At=0.5 s, Ax=2 
cm y T=120 s, la precisión de la solución numérica reportada en (Millán et al., 2011) 
presenta un error absoluto acumulado de L1=3.64 respecto a la solución analítica reportada 
en el mismo artículo. 

A continuación se presentan las gráficas de: las condiciones iniciales, la solución numérica 
reportada en este artículo y la solución analítica reportada. 


Ciencia e Ingeniería 2014 Vol. (1) N* (1) Lyon E. y Rosales L. 
FUNDAMENTOS DE SIMULACIÓN MEDIANTE ALGORITMO EN DIFERENCIAS FINITAS 
DE LA EVOLUCIÓN ESPACIAL Y TEMPORAL EN LA CONDUCCIÓN DEL CALOR 


GRÁFICA 1: SOLUCIÓN NUMÉRICA REPORTADA EN [1] 


100 y 7 


Solución numérica reportada en [1] Dt=0.5s Dx=2cm T=120s 
r r A 7 7 


Patrón inicial 
de temperatura [1] | 


7o0+ a Solución >< z 
A numérica[1] 


30+ Y Solución ? Ñ a 
Analítica [1] Ñ 


ñ 
0) 2 4 6 5] 10 12 14 16 18 20 


En la GRÁFICA 2 puede observarse que para las condiciones impuestas de At=0.5 s, Ax=2 
cm y T=120 s, la precisión del algoritmo reportado en este artículo presenta un error 
absoluto acumulado de L1=3.69 respecto a la solución analítica reportada en (Millán et al., 
2011). 

Comparando los errores acumulados se observa una diferencia de 1.37 %, entre los 
resultados reportados en (Millán et al., 2011) y los reportados en el presente artículo, 
tomando los primeros como patrón de referencia. No obstante de no conocer los criterios de 
redondeo al presentar los valores en (Millán et al., 2011) podemos afirmar que la diferencia 
entre los valores de L1 es aceptable. 


GRÁFICA 2: SOLUCIÓN NUMÉRICA REPORTADA EN ESTE ARTÍCULO 


Solución numérica para validar artículo Dt=0.5s Dx=2cm T=120s 
100 Y Tr 7 y Y 


A r Y - 
/ 


A Patrón inicial 
sob PS Ñ de temperatura [1] | 


7oL 4 Solución y) 
numérica Artículo 


Solución 
Analítica [1] 


10| Y LA = 3.69 SN 
9d 1 m7 1 1 mi L 1 L 


1 
o 2 + 6 8 10 12 14 16 18 20 
xcem 


3. Afinación de la ventana temporal de simulación (T). 


Se estudió el comportamiento del algoritmo al afinar la ventana temporal de simulación 
(T) dejando fijos At y At, midiendo el acumulado L1. 


A continuación se muestran las gráficas con los resultados de afinación de la ventana 
temporal de simulación T para seis valores 25, 50, 75, 100, 125 y 150, se midió L1 para 
cada caso. 


En la GRÁFICA 3 se presentan los resultados de la solución numérica para seis valores 


de la ventana temporal de simulación T, analizando los valores de T y L1 podemos afirmar 
que el óptimo se encuentra entre las ventanas T=100 y T=150 muy cercano a 125, la 


20 
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solución numérica para este valor en la gráfica casi se solapa con la solución analítica. No 
se presentó el fenómeno de la inestabilidad. 


GRÁFICA 3: AFINACIÓN DE LA VENTANA TEMPORAL DE 


SIMULACRO (T) 


Afinación de la ventana de tiempo T con Dt=0.5 s Dx=2 cm 
100 " - - a E ; 


Terperaua”c 
J 


4. Afinación del paso Ax. 


Se estudió el comportamiento del algoritmo al afinar del paso Ax dejando fijos At=0.5 y 
T=125. 


La prueba se diseñó tomando en cuenta la condición de estabilidad enunciada en 
(Ezquerro, 2012) 0 <A < 0.5 (p. 98). En la GRÁFICA 4 se muestran los resultados de 
afinación del paso en x para cinco valores Ax [0.55, 0.50, 0.45, 0.40] menores que 1 y 
aproximarnos al límite de estabilidad 1=0.5, A [0.25, 0.30, 0.37, 0.47] para observar su 
influencia en la convergencia y en la estabilidad. 


En la GRÁFICA 4 se observa que las soluciones numéricas para los cinco valores Ax < 
1 se solapan, lo cual indica que la influencia sobre el error es despreciable, esto es debido a 
que en este algoritmo el error por truncamiento depende de (Ax)”, no aparece el fenómeno 
de la inestabilidad pues no hemos sobrepasado el límite 1=0.5 


GRÁFICA 4: PASO Ax PARA COMPORTAMIENTO ESTABLE 2<0.5 


Afinación del paso Dx con T=125 y Dx=0.5 Y<0.5 
100 Y 7 7 Y Y Y 


TemperdurarC 
h 16) 
o [2] 
ES 
1 ] 


Y > 
30( Pa Dx Y Sa e 4 
A 0.55 0.25 
20F 9 0.50 0.30 DAS 7 
á 0.45 0.37 Xx 
10t Eh 0.40 0.47 SN 
es 4 Ñ 
(o) 2 4 6 8 10 12 14 16 18 20 
xcm 


En la GRÁFICA 5 se observa que la solución numérica para los valor Ax =0.35 y A= 
0.65, entra en inestabilidad tal y como lo plantea lo enunciado en (Ezquerro, 2013) para la 
condición de estabilidad O < A < 0.5 (p. 98) 
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GRÁFICA 5: PASO Ax PARA COMPORTAMIENTO INESTABLE 2<0.5 


x 10% Afinación del paso Dx con T=125 y Dx=0.5 Y>0.5 
8 7 Y Y Y y Y Y y 


Terpsaua rc 
o 


5. Afinación del paso At. 


Se estudió el comportamiento del algoritmo al afinar del paso At dejando fijos Ax=2 y 
T=125. 


La prueba se diseñó tomando en cuenta la condición de estabilidad enunciada en 
(Ezquerro, 2012) 0 < 1 < 0.5 (p. 98). En la GRÁFICA 6 se muestran los resultados de 
afinación del paso en t para tres valores At [12, 13, 14] y aproximarnos al límite de 
estabilidad 1=0.5, 1 [0.45, 0.49, 0.53] para observar su influencia en la estabilidad. 


En la GRÁFICA 6 se observa que la solución numérica para los valores At = 13 y A= 
0.49 comienza a entrar en inestabilidad y para At = 14 y A= 0.53, ya el fenómeno de 
estabilidad es notable como lo plantea lo enunciado en (Ezquerro,2012) para la condición de 
estabilidad O < 4 < 0.5 (p. 98). Se puede notar que el algoritmo es menos sensible al 
acercarnos al límite de estabilidad cuando variamos At, comparado con el caso de la 
afinación de Ax. 


GRÁFICA 6: PASO At ALREDEDOR DEL LÍMITE DE ESTABILIDAD 1-0.5 


Afinación del paso Dt con Dx=2 T=125 
100 r 7 eN r 


Temperatura ”C 
E] 
o 0 


En la GRÁFICA 7 se observa el comportamiento temporal para los nodos en x y se 
observa como ser propaga la inestabilidad en el tiempo, este es un fenómeno que no es 
característico del sistema en estudio. Se da bajo simulación cuando la razón de muestreo 


22 
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AU(Ax)? no es apropiada respecto al coeficiente de difusión térmica a, lo que termina 
violando la condición de estabilidad del algoritmo 0<1<0.5. 


GRÁFICA 7: COMPORTAMIENTO TEMPORAL DE LA INESTABILIDAD EN 


LOS NODOS 
Comportamiento temporal de la temperatura en los nodos lado izquierdo 
so — - 
70H llo A 7 
so 33 a =l 
AAA A o A E Nodo 5 
Ó so A a A y 
rd a 
— == 
40 ds E 
¿ o E Nodo 4 
E 30| PS ] 
Nodo 3 
20 E E 7 
Nodo 2 
10H A 
Nodo 1 
o a 
o 20 40 60 80 100 120 
ts 
VIH. CONCLUSIONES. 


El método explícito basado en la aproximación de la derivada temporal mediante diferencia 
progresiva y la derivada espacial mediante diferencia de punto central, es fácil de implementar, 
pero requiere de un conocimiento suficiente de la dinámica del sistema a simular para evitar 
llevar el algoritmo a la inestabilidad por un muestreo o paso de simulación inadecuado. 


El modelado mediante ecuaciones diferenciales parciales parabólicas permite estudiar el 
fenómeno en su evolución temporal y espacial. 


La dinámica del sistema determinada por las características de difusión térmica del material, 
conjuntamente con la razón de paso temporal espacial configura la base para el análisis de 
convergencia y estabilidad del algoritmo. 


El algoritmo planteado en este artículo se comporta aceptablemente al compararlo con los 
resultados del artículo base (Millán et al., 2011). 


Para el caso de estudio la convergencia del algoritmo varía con la ventana temporal. 
Para el caso de estudio la estabilidad del algoritmo no es sensible a la variación de la ventana 
temporal. 


Durante el refinamiento de los pasos temporal y espacial la estabilidad el algoritmo mostró 
mayor sensibilidad cuando variamos el paso espacial, aunque durante la variación del paso 
temporal también se presenta la inestabilidad pero con menor sensibilidad. 
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